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ABSTRACT 

As  part  of  the  international  Greenland  Sea  Project,  Woods  Hole  Oceanographic 
Institute  and  Scripps  Institution  of  Oceanography  deployed  a  six  transceiver  ocean 
acoustic  tomography  array  to  monitor  ocean  ventilation  and  circulation  over  the 
19SS-89  winter  cooling  season.  A  stochastic  inverse  method  computer  code  which  at- 
tains a  solution  by  minimizing  mean  square  error  is  used  to  perform  inversions  of  the 
Greenland  Sea  tomography  data. 

A  computer  simulated  ocean  is  used  to  evaluate  various  aspects  of  system  perform- 
ance. We  first  consider  the  advantages  and  problems  associated  with  using  a  ray  theory 
based  algorithm  to  establish  the  forward  problem  for  the  Greenland  Sea  tomography 
array.  Next,  we  made  two  adjustments  to  our  inversion  code  and  discuss  the  effects  on 
system  performance.  The  first  adjustment  allows  for  layers  of  different  thicknesses  in  the 
inverse  solution  to  increase  the  density  of  estimates  in  regions  of  interest.  The  second 
adjustment  allows  the  estimator  to  expect  variability  of  the  unknown  field  to  decrease 
exponentially  with  depth. 

Our  results  show  the  ray  theory  based  algorithm  is  an  adequate  method  of  modeling 
ray  paths  in  the  Greenland  Sea,  but  has  limitations.  Reliability  of  ray  paths  degrades 
as  launch  angles  become  shallower  and  if  strong  gradients  and  rapidly  changing  gradi- 
ents in  sound  speed  are  present  in  the  vicinity  of  the  transceiver  elements.  We  also  find 
varying  the  thickness  of  layers  in  the  solution  allows  us  to  examine  the  more  variable 
upper  ocean  in  greater  detail  without  increasing  computational  effort.  However,  this 
adjustment  alone  has  the  undesirable  side  effect  of  shifting  system  resolution  towards  the 
lower  ocean.  This  shift  in  resolution  is  offset  by  informing  the  estimator  about  the  ver- 
tical variability  distribution  of  the  unknown  field.  This  a  priori  knowledge  is 
parameterized  by  the  covariance  function  of  the  unknown  field.  Uncertainty  in  knowing 
the  true  variability  distribution  affects  model  performance.  The  inverse  solution  is  more 
sensitive  to  underestimating  than  overestimating  the  true  value  of  folding  depth.  The 
model  is  also  more  sensitive  to  both  underestimation  and  overestimation  at  small  true 
folding  depths. 

A  set  of  Greenland  Sea  data  between  one  tranceiver  pair  was  processed  by  Woods 
Hole  Oceanographic  Institute.  Although  only  three  groups  of  eigenrays  are  involved, 
initial  inversion  results  indicate  the  estimator  detects  seasonal  changes  and  synoptic  scale 
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events  occurring  at  time  scales  greater  than  20  days,  however,  solutions  show  wide  fluc- 
tuations at  time  scales  shorter  than  20  davs. 
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I.     INTRODUCTION 

A.     OCEAN  ACOUSTIC  TOMOGRAPHY 

Tomography  is  a  method  to  remotely  sense  interior  structures.  It  uses 
electromagnetic,  seismic,  or  sound  waves  to  propagate  through  and  probe  relatively 
transparent  media.  For  example,  x-rays  are  commonly  used  in  medical  computer  as- 
sisted tomography  (CAT)  scans  and  seismic  waves  induced  by  surface  sources  are  used 
to  study  the  interior  of  the  earth  (Backus  and  Gilbert, 1967).  Using  analagous  methods, 
Munk  and  Wunsch  (1979)  introduced  ocean  acoustic  tomography  as  a  means  of  moni- 
toring mesoscale  fluctuations  in  ocean  basins.  Sound  waves,  which  are  sensitive  to 
density  changes  and  currents  but  transit  well  through  the  ocean,  gather  information  as 
they  encounter  structures  in  the  ocean  interior.  The  data  is  received  in  the  form  of  a 
peturbated  travel  time  from  a  source  to  a  receiver.  Using  inverse  techniques,  the  best 
estimate  as  to  what  structure  could  result  in  the  observed  data  is  constructed. 

Ocean  acoustic  tomography  has  several  practical  advantages  over  conventional 
hydrographic  study  techniques,  as  cited  by  Chiu  ei  al.  (19S7).  The  monitoring  system 
can  be  established  as  a  semipermanent,  continuous  observing  system  which  is  not  greatly 
affected  by  weather.  It  has  high  temporal  resolution  and  can  cover  an  extensive  volume 
of  ocean  simultaneously  with  only  a  few  moorings,  which  reduces  deployment  and 
maintenance  costs.  Another  advantage  stated  by  Munk  and  Wunsch  (1979)  is  that  the 
amount  of  additional  information  gained  by  each  additional  mooring  increases  more 
rapidly  than  that  gained  by  additional  "spot"  measurements  since  each  mooring  sets  up 
distinct  new  paths  to  each  of  the  previously  set  moorings.  The  gain  in  spot  measure- 
ments is  more  or  less  linear  with  the  number  of  instruments  deployed. 

Ocean  acoustic  tomography  has  been  successfully  demonstrated  in  a  variety  of  ap- 
plications since  its  introduction.  Some  examples  include  mapping  mesocale  midocean 
eddy  fields  (Comuelle  el  al..  19S5).  analysis  of  planetary  waves  (Chiu  ei  «/.,19S7), 
studying  of  ocean  currents  (DeFerrari  ei  al.,  1986),  and  surface  wave  analysis  (Lynch  el 
a/..19S7).  It  has  also  been  chosen  as  a  monitoring  tool  as  part  of  the  Greenland  Sea 
Project.    Its  application  in  the  Greenland  Sea  is  the  impetus  for  this  thesis. 

The  Greenland  Sea  is  a  well  known  contributor  to  world  ocean  ventilation.  Several 
conceptual  theories  have  been  developed  to  describe  the  process,  but  direct  observations 
of  ventilation  are  scarce.    The  Greenland  Sea  Project  was  put  forward  by  the  Arctic 


Ocean  Sciences  Board  to  gain  understanding  of  the  processes  relevant  to  deep  con- 
vection (Meincke.1989). 

B.  THE  GREENLAND  SEA  PROJECT 

The  Greenland  Sea  Project  (GSP)  is  a  five  year  (1987-1992).  international  scientific 
study  to  understand  the  large  scale,  long  term  interactions  of  air,  sea,  and  ice  in  the 
Greenland  Sea.  The  five  key  elements  of  the  GSP  are:  (1)  a  study  of  air-sea-ice  inter- 
action. (2)  a  study  of  ocean  ventilation.  (3)  a  study  of  ocean  circulation  and  mixing,  (4) 
a  study  of  atmospheric  energetics,  and  (5)  a  study  of  biological  processes.  The  role  of 
ocean  acoustic  tomography  is  to  assist  in  the  studies  of  ocean  ventilation  and  circu- 
lation. The  array  of  transceiver  moorings  is  designed  to  measure  changes  in  the  inte- 
grated properties  of  the  Greenland  Sea  gyre  through  a  winter  cooling  season  and  may 
also  provide  valuable  information  in  measuring  the  response  of  the  wind  driven  gyre  to 
changes  in  the  curl  of  the  wind  stress  (Greenland  Sea  Science  Planning  Group.  1986). 

A  six  transceiver  array  was  successfully  deployed  in  the  Greenland  Sea  during 
September-October.  19S8.  Figure  1  depicts  the  planned  deployment  area  and  geometry 
of  the  array.  The  moorings  were  shifted  slightly  from  their  original  planned  location  to 
accommodate  for  the  rough  bathymetry  found  at  the  original  Mooring  3  site.  Mooring 
2  was  redeployed  due  to  leakage  of  some  O-ring  seals.  The  redeployed  mooring  is  des- 
ignated here  as  Mooring  2a.  The  geodic  coordinates  of  the  moorings,  the  source  depths, 
and  the  receiver  depths  are  given  in  Table  1.  (Worchester  and  Howe,  1989). 

C.  THESIS  OBJECTIVES 

There  are  three  basic  objectives  of  this  thesis.  The  first  is  to  discuss  the  advantages 
and  problems  of  predicting  ray  paths  between  array  transceivers  using  an  algoritm  based 
on  ray  theory.  The  Greenland  Sea  offers  an  acoustic  environment  which,  in  conjunction 
with  the  array  geometry,  makes  ray  path  determination  challenging  to  model.  Accurate 
ray  path  prediction  is  an  essential  element  of  establishing  the  "forward"  problem. 
Chapter  III  contains  our  discussion  on  this  topic. 

The  second  objective  is  to  develop  an  inversion  code  appropriate  for  application  to 
Greenland  Sea  tomography  data.  We  begin  with  an  inversion  code  originated  by  Chiu 
and  modified  by  Kao  (1989).  The  code,  or  estimator,  treats  the  ocean  as  a  volume  which 
has  been  subdivided  into  discrete  boxes.  The  boxes  have  equal  horizontal  dimensions 
and  the  vertical  layers  of  the  ocean  are  equally  spaced.  An  estimate  of  the  unknown 
variable  is  calculated  for  each  box.   The  estimator  also  treats  the  ocean  as  a  statistically 


Figure   1.       Planned  Deployment  Site  of  (he  Greenland  Sea  Tomography  Array 


Table   1.     TRANSCEIVER  LOCATIONS 


Mooring 

Lat-Long  Position 

X-Y  Position  (km) 

Source 
Depth  (m 

Receiver 
Depth  (in 

1 

75°58.08'A\(X)r50.00'Jr 

154.790   220.317 

99.7 

150.4 

2 

75°03.69'Ar,000°40.25'£ 

224.683    119.080 

94.0 

144.7 

2a 

75°03.88'A\000°38.20'£ 

223.715    119.491 

97.4 

148.1 

3 

74o09.38'A".001°52.90'W' 

148.972     18.146 

96.0 

146.7 

4 

74°28.90'A\005°47.30'  W 

33.270    60.958 

94.6 

117.7 

5 

75°34.27'A\006°07.70'ir 

34.547    1S2.902 

101.5 

124.6 

6 

75e03.60'A'.002°58.00'n' 

120.000    120.000 

92.5 

143.2 

homogeneous  volume.  Therefore  it  expects  variability  (RMS  values  of  the  unknown 
variable)  in  any  region  of  the  ocean  to  be  the  same. 

We  have  made  two  fundamental  adjustments  to  the  inversion  code  and  discuss  how 
they  effect  system  performance.  Our  first  adjustment  is  to  vary  the  layer  thickness  so 
that  there  is  a  greater  density  of  boxes  in  the  region  of  interest.  This  adjustment  allows 
for  more  detail  in  the  region  of  interest  without  increasing  the  computational  effort. 
Our  second  adjustment  tells  the  estimator  to  expect  a  depth  dependent  variability  dis- 
tribution of  the  unknown  variable.  In  this  way.  the  estimator  produces  solutions  which 
are  more  realistic  statistically.  The  effects  of  the  adjustments  on  system  performance  is 
evaluated  by  conducting  inversions  on  computer  simulated  oceans  with  known  statistical 
parameters.   The  results  of  these  studies  are  also  contained  in  Chapter  III. 

For  convenience  in  coding,  we  use  a  Cartesian  coordinate  system.  Based  on  ranges 
computed  for  the  WGSS4  spheroid  (Worchester  and  Howe,  1989).  we  converted  the 
latitude  and  longitude  positions  to  an  XV  plane  which  closely  approximates  the  geodic 
orientation.  The  plane  is  240  km  by  240  km  and  centered  on  Mooring  6.  The  XY  co- 
ordinates are  included  in  Table  1.  Though  there  is  some  slight  distortion  in  the  overall 
shape  of  the  array,  the  error  between  the  geodic  ranges  and  the  model  ranges  is  less  than 
15  meters  along  each  path,  as  shown  in  Table  2  below. 


Table  2.  MOORING  RANGES  (KM)  AND  RANGE  ERRORS  (M):  Ranges 
(above  diagonal)  are  based  on  the  WGSS4  spheroid,  range  errors  (below 
diagonal)  are  XV  ranges  minus  spheroid  ranges. 


Mooring 

1 

2 

2a 

3 

4 

5 

6 

1 

0 

123.033 

122.120 

202.258 

200.405 

125.922 

106.179 

2 

-13 

0 

1.039 

126.177 

200.057 

200.576 

104.701 

2a 

13 

13 

0 

125.927 

199.232 

199.500 

103.706 

3 

-3 

-3 

-1 

0 

123.367 

200.594 

105.896 

4 

1 

-14 

5 

2 

0 

121.947 

104.921 

5 

8 

-14 

13 

-1 

4 

0 

106.109 

6 

-1 

-14 

10 

_2 

_2 

-1 

0 

Our  third  objective  is  to  test  the  inversion  code  with  observed  data  derived  from  the 
Greenland  Sea  tomographic  array.  The  studies  on  system  performance  discussed  earlier 
are  conducted  using  computer  simulated  oceans  with  known  characteristics.  However, 
the  simulated  oceans  are  limited  in  the  scope  to  which  they  can  minimi ck  the  actual 
conditions  found  in  the  Greenland  Sea.  Using  actual  Greenland  Sea  travel  time  data  is 
the  only  way  to  determine  if  the  model  response  will  be  realistic.  Currently  there  is  a  set 
of  data  derived  from  signals  transmitted  between  Mooring  4  and  Mooring  5  available  for 
input  into  our  model.  The  data  set  was  processed  by  Woods  Hole  Oceanographic  In- 
stitute (YYHOI).  Though  this  will  not  enable  us  to  form  significant  conclusions  about 
the  array  site  from  such  a  small  data  base,  ingesting  real  data  into  our  model  is  helpful 
in  preparing  the  inversion  code  for  the  complete  data  set  when  it  becomes  available. 
Our  results  are  discussed  in  Chapter  IV. 

To  complete  our  study,  we  provide  the  following  information.  In  Chapter  II,  we 
present  an  overview  of  the  basis  on  which  our  system  is  developed.  Here,  the  math- 
ematical foundation  of  the  tomographic  "forward"  and  "inverse"  problems  is  discussed. 
Then  in  Chapter  V.  we  summarize  our  final  conclusions  and  make  recommendations  for 
future  improvements  related  to  this  work. 


II.     OVERVIEW  OF  ACOUSTIC  TOMOGRAPHY  PRINCIPLES 

In  general,  the  application  of  tomogrgaphy  to  monitoring  the  ocean  interior  is  done 
in  two  parts.  In  the  "forward"  problem,  we  establish  a  model  describing  the  physical 
relationship  between  the  observed  data  and  the  unknown  variables  which  are  to  be  esti- 
mated. We  will  cast  the  forward  problem  as  a  Fredholm  Integral  Equation  of  the  First 
Kind  in  the  form 

4=  I      Si(s)As)ds  +  et        i-  l,2,...,m  (2.1) 

''path 

where  d,  is  the  data  accumulated  along  the  /'*  path,  #,(5)  is  the  data  kernel  which  contains 
the  physics  and  operates  on  the  unknown,  f[s)  is  the  unknown  structure  which  is  to  be 
estimated,  and  e,  is  the  experimental  noise  accumulated  along  the  /"•  path. 

The  "inverse"  problem  is  to  determine  a  solution  to  the  unknown  structure  which  fits 
the  data.  An  undetermined  problem,  one  generally  finds  many  solutions  that  fit  the  data 
and  must  choose  an  optimal  solution  based  on  some  objective  criteria.  For  the 
Greenland  Sea  transceiver  array,  the  observed  data  are  travel  times  of  sound  energy 
along  reciprocal  acoustic  multipaths  which  "connect"  sources  and  receivers.  From  these 
data,  we  can  estimate  both  sound  speed  perturbations  with  respect  to  a  reference  ocean 
(density  tomography),  and  ocean  currents  (current  tomography). 

A.     THE  FORWARD  PROBLEM 

Sound  propagation  in  the  ocean  is  a  function  of  temperature,  pressure,  and  salinity 
as  well  as  a  function  of  the  current  field.  For  the  forward  problem,  we  need  to  determine 
the  acoustic  multipaths  which  sound  energy  follows  from  source  to  receiver.  To  model 
these  "eigenrays",  we  use  an  algorithm  based  on  ray  theory.  Ray  theory  provides  a 
simple  physical  description  of  the  acoustic  multipaths  and  the  modeling  equations  are 
straightforward.   We  start  with  the  eikonal  equation, 

I  VO(W)  |2  =  n(xytz?  (2.2) 

where  O  is  the  acoustic  phase  defining  a  wave  front  and  n  is  the  index  of  refraction. 
Following  the  example  of  Ziomek  (1985),  we  neglect  horizontal  refraction  so  that  the 
predicted  ray  paths  lie  on  a  vertical  plane  passing  through  the  source  and  receiver. 


Following  the  derivation  used  by  Kao  (1989).  and  treating  the  reference  ocean  as  a 
motionless  medium,  we  arrive  at  the  key  equation  for  modeling  ray  paths  in  the 
Greenland  Sea. 


(2.3) 


where  z  is  the  depth  of  the  ray  at  range  R,  6  is  the  ray  path  angle  with  respect  to  the 
horizontal,  and  the  subscript  o  indicates  conditions  at  the  acoustic  source. 

The  speed  of  the  wave  front  through  the  true  ocean  is  affected  by  both  the  sound 
speed  and  current  structure  along  the  path.   We  express  this  as 

:   =c{xy,z,t)  +  V{xj,z,t).s  (2.4) 


dl 


where  c  is  sound  speed.  V  is  the  ocean  current  field,  and  5  is  a  unit  vector  tangent  to  the 
ray  path.  5  .  and  pointed  in  the  direction  of  sound  energy  propagation.  From  Equation 
(2.4),  we  can  determine  the  sound  travel  time  along  an  eigenray  path. 


T  = 


Jpath\   c  +  V  •  s    / 


We  further  separate  c  into  the  sum  of  sound  speed  in  a  reference  ocean  and  sound  speed 
perturbations  with  respect  to  this  reference, 

c  =  c0  +  be  (2.6) 

Travel  time  can  also  be  separated  into  the  sum  of  travel  time  along  an  eigenray  in 
a  reference  ocean  and  the  perturbation  in  travel  time  due  to  effects  in  the  ocean  not  re- 
presented by  the  reference  conditions. 

T=T0  +  ST  (2.7) 

Substituting  Equations  (2.6)  and  (2.7)  into  Equation  (2.5),  we  have  '  '■ 


T0  +  dT  = 


path 


1  + 


Sc  +  V 


ds 


(2.8) 


Under  the  mild  restriction  that  |  dc  4-  V  •  5  |  <  <  |  cc  |  .we  can  neglect  the  higher  order 
terms  of  a  Taylor  series  expansion  of  the  term  in  parenthesis  and  arrive  at 


T0  +  STc 


path 


iO-^v. 


(2.9) 


In  density  tomography,  which  uses  one  way  travel  times,  the  effect  of  ocean  currents 
is    much    smaller    than    the    elTects    of    temperature,    pressure,    and    salinity,    i.e. 


I  V  •  5  I  <  <  I  Sc  I  .    Since  we  have  defined  T„  =        ~  ds,  Equation  (2.9)  reduces  to  the 

path     ° 

linearized  form 


ST 


*  path 


-  1 


[Sc)ds 


(2.10) 


In  current  tomography,  we  must  consider  reciprocal  paths.  The  effects  of  temper- 
ature, pressure,  and  salinity  are  the  same  in  either  direction  whereas  the  effects  of  the 
current  field  will  be  opposite  under  these  circumstances.  Therefore,  imbedded  informa- 
tion about  the  currents  may  be  extracted  by  subtracting  the  two  reciprocal  path 
equations. 

Consider  the  the  travel  time  for  sound  along  a  path  indicated  by  =>  calculated  from 
Equation  (2.9). 


t0  +  st- 


path 


-Ml- 


dc  +  V  •  S 


ds 


(2.11) 


and  for  its  reciprocal  path  indicated  by 


T0  +  dT 


Cn 


Sc  +  V  •  £ 

Co 


pam 


:/.9 


(2.12) 


Since  V  •£-  =  —  ¥•£**  at  all  points  on  the  paths,  subtracting  Equation  (2.12)  from 
Equation  (2.1 1 )  yields  a  linearized  expression  for  current, 


3T=  -(sr 


-  1 


;v .  rys 


(2.13) 


'path 


Thus  far.  the  forward  model  has  treated  travel  time  data  only  as  noise  free  meas- 
urements. In  practice,  experimental  noise  corrupts  the  signal  and  contributes  to  the  er- 
ror in  the  estimates.  Eactors  such  as  tidal  forces  which  cause  displacement  of  transceiver 
elements  and  internal  waves  are  sources  of  experimental  noise.  Signal  processing  is  used 
to  minimize  these  effects.  Details  of  signal  processing  techniques  used  for  ocean 
tomography  can  be  found  in  Spindel  (1985). 

Including  experimental  noise  in  our  model,  we  now  have  expressions  which  relate 
the  unknown  variables  to  practical,  observable  data.    Eor  density  tomography,  we  have 


ST= 


path 


-1 


ldc)ds  +  e6c 


(2.14) 


and  for  current  tomography,  we  have 


sT~-sr 


••path 


(V  .  s)ds  +  c- 


(2.15) 


The  similarity  of  the  data  kernels,  which  contain  the  physics  of  the  model,  shows  the 
close  relationship  between  the  two  areas  of  tomography. 


B.     THE  INVERSE  PROBLEM 

Once  the  forward  problem  has  been  established,  we  search  for  unknown  structures 
which  fit  the  observed  data.  Since  we  cannot  sufficiently  sample  the  ocean  to  have  a 
complete  set  of  independent  data,  the  problem  is  underdetermined,  with  an  infinite 
number  of  possible  solutions  that  will  fit  the  data.  To  reach  a  unique  solution,  we  must 
devise  an  estimator  which  imposes  additional  constraints  on  the  system.  Choosing  an 
estimator  based  on  the  Gauss-Markoff  Theorem  allows  the  use  of  statistical  information 
to  evaluate  system  performance. 
1.     The  Estimator 

"We  approximate  the  continuous  system  in  a  discretized  form.  This  allows  us  to 
cast  the  problem  into  matrix  algebra  for  straightforward  implementation  of  the  estimator 
in  a  numerical  computer.  Applying  the  Gauss-MarkofT  Theorem  (Liebelt,1967)  and 
following  the  works  of  Cornuelle  a  al.  (1985)  and  Chiu  et  al.  (1987),  our  criterion  for  the 
"best"  solution  of  the  unknown  structure  is  one  that  is  linear  with  the  data  and  has 
minimal  mean  square  error  with  respect  to  the  true  solution.  Additional  contraints  on 
the  system  are  imposed  through  the  specification  of  an  a  priori  covariance  matrix. 
Posing  the  forward  problem  in  discretized  form,  we  have 

d=Af+e  (2.16) 

where  the  m  dimensional  column  vector,  d,  contains  the  data  derived  from  travel  times 
of  m  resolvable  eigenrays.  The  n  dimensional  column  vector,/",  is  the  unknown  structure 
to  be  estimated.  The  ocean  has  been  divided  into  n  discrete  boxes,  within  each  the  un- 
known variable,  3c  or  V.  is  assumed  to  be  constant.  In  our  treatment,  the  discrete  boxes 
have  equal  dimensions  in  the  horizontal  plane,  forming  squares.  However,  the  vertical 
dimension  varies  with  depth,  forming  uniform  layers  of  different  thicknesses.  The  ex- 
perimental noise,  e .  which  contaminates  the  data  along  each  ray  path,  is  also  m  dimen- 
sional. The  linear  operator  matrix,  A,  has  m  x  n  dimensions.  The  rows  of  the  operator 
matrix  are  analagous  to  a  set  of  data  kernels  of  the  continuous  case.  They  contain  the 
physics  which  relates  the  data  and  unknown  fields. 

From  the  Gauss-MarkofT  Theorem,  the  optimal  estimate  of  an  unknown  pa- 
rameter can  be  calculated  from 

f  =  Ce\TC;]J  (2.17) 

where  the  covariance  matrix  of  the  total  error,  c.  of  the  estimate  is 
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Cf  =  Cf  -  (C{AT)(ACjAT+  C,r,(C/Ar)r  (2.18) 

The  diagonal  elements  of  C,  are  the  minimum  mean  square  errors  of  the  estimate.  C,  is 
the  covariance  matrix  of  experimental  noise,  which  we  assume  is  known,  and  C,  is  the 
covariance  matrix  of  the  unknown  variable.  It  is  through  Cf  that  we  insert  a  priori  in- 
formation in  the  estimator.  To  reach  an  optimal  solution,  we  assume  the  unknown  fields 
have  Gaussian  distributions  and  are  statistically  homogeneous  throughout  each  layer. 
We  express  these  conditions  through  a  covariance  function,  Cf,  which  is  cast  into  matrix 
form  as  C, .   The  covariance  function  is  expressed  as, 


Cf=  o(x',y',  z')a{x",y",  2")  exp 


I 


V      / 


*-"*  x      J  \        ^~*\*       I  *       ^^  ^ 


(2.19) 


where  the  separation  of  two  points  in  each  dimension  is  given  as 
A.v  =  x"  —  x',  Ay  —y"  —  y' ,  and  Az  =  2"  —  2',  respectively.  The  dimensional  correlational 
lengths.  L,.  Lv.  and  L,  tell  the  estimator  how  well  information  at  one  point  relates  to 
surrounding  points  in  the  system.  It  is  optimal  to  use  the  true  correlation  lengths,  but 
the  true  values  may  not  be  well  known.  To  avoid  overconstraining  the  system,  i.e.  to 
avoid  correlating  information  which  is  not  truly  correlated,  we  can  use  conservatively 
small  values  of  Lx.  Ly,  and  L..  We  characterize  the  vertical  variability  distribution  of  the 
unknown  field  through  a.  which  is  treated  by  the  estimator  as  an  exponentially  decaying 
function  with  depth,  parameterized  by  a  surface  value  and  an  e  -folding  depth.  In  each 
layer,  an  average  value  based  on  this  curve  is  used.  With  a  priori  knowledge  in  place, 
the  inversion  code  can  generate  a  unique  least  mean  square  error  solution  from  observed 
data  as  well  as  error  bars  and  system  resolution  measures. 
2.     Measurement  of  Error 

The  total  error  in  the  estimate  can  be  cast  as  two  statistically  independent  terms; 
bias,  b  =  <f>  -f.  and  random  error,  Af=f—  <f>  .  The  covariance  of  error  is  simply 
the  sum  of  its  components, 

C£  =  <Mr>  +  C^  (2.20) 


where 


CA;  =  CE  A^J1  AC£  ■  (2.21) 
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Following  Chiu  et  al.  (19S7)  and  Wiggins  (1972).  we  apply  a  powerful  tool  called  singular 

value  decomposition  of  matrices.   The  linear  operator  matrix.  A,  is  scaled  and  factorized 

into  a  diagonal  matrix  of  eigenvalues,  A,  and  two  matrices  of  associated  eigenvectors. 

U  and  V  such  that. 

_j_       j_ 
C£    2AC/  =  UAVr  (2.22) 

Each  diagonal  element  of  A.  /,,  has  two  associated  eigenvectors,  u,  and  v, ,  which  are  the 
i"'  columns  of  U  and  V,  respectively.     The  eigenvectors  within  U  and  V  are  sets  of 
orthogonal  basis  vectors  in  the  data  space  and  unknown  parameter  space,  respectively. 
Applying  matrix  decomposition  as  in  Chiu  et  al.  (1987),  we  now  write  the 
covariances  of  total  error  and  random  error  as 

_L  _L 

C£  =  C/[I  -  VA(I  +  A2)_1AVr]  C/  (2.23) 

and 


j 

,  .  -l  =  V-y   L  ▼  \*  T  /V  )      J\   {l  T  J\    I       »      J  K,f' 


C  ;  =  C,-2[V(I  +  A2r'A2(I  +  A2)~1\'7"]  Qf2  (2.24) 


where  I  is  an  n  x  n  identity  matrix.  We  derive  an  expression  for  bias  by  taking  the  ex- 
treme case  where  all  eigenvalues  approach  infinity  so  that  Q^  aproaches  zero.  From 
Equation  (2.23)  we  have, 

j_  j_ 

<bhT>  =  C/(I  -  Wr)C/.  (2.25) 

3.     Measurement  of  Resolution 

Besides  a  measure  of  error,  we  also  evaluate  system  performance  through  its 
ability  to  resolve  ocean  features.  Continuing  with  singular  value  decomposition  of  ma- 
trices, we  obtain  an  expression  for  the  optimal  estimate  in  terms  of  eigenvalues, 


(M/Ce    V)v;  (2.26) 


where  k  is  the  number  of  non-zero  eigenvalues.   Equation  (2.26)  shows  the  solution  is  a 
weighted  sum  of  eigenvectors  where  the  eigenvalues  control  the  weighting.     The  ex- 
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pression  in  parenthesis  represents  the  components  of  the  data  vector  which  has  been 
expanded  into  a  linear  combination  of  weighted  eigenvectors.  As  demonstrated  by  Chiu 
ei  al.  (19S7).  the  )}x  values  are  analagous  to  signal-to-noise  ratios  for  each  component 
of  the  expanded  data  vector.  For  /?  >  >  1,  signal  dominates  noise  in  that  particular 
component,  and  the  reverse  is  true  for  X]  <<  1.  From  Equation  (2.26)  we  see  that  the 
effect  of  noise  on  the  estimate  is  minimized  by  downweighting  those  components  asso- 
ciated with  high  noise  levels. 

We  now  define  the  resolution  matrix,  R,  which  has  n  x  n  dimensions, 

R  =  VA(I  +  A2)_1AVr  (2.27) 

The  i'h  row  of  R  is  the  resolution  kernel  of  the  /'*  box  in  the  discretized  ocean.  Substi- 
tuting into  Equation  (2.23).  we  get 

_L  JL 

CE  =  C/(I  -  R)C/  (2.28) 

Ideally,  if  the  resolution  matrix  is  an  identity  matrix,  the  system  has  no  error.  In  the 
limit,  as  we  approach  the  continuous  case  with  complete  and  noise  free  data,  the  resol- 
ution kernels  approach  delta  functions  and  the  system  would  be  perfectly  resolved.  In 
practice,  however,  the  resolution  kernels  have  side  lobes  and  amplitudes  less  than  unity 
along  the  diagonal.  In  Equation  (2.28),  the  relationship  between  resolution  and  error  is 
linear.  Therefore,  we  may  use  the  resolution  kernel  peak  value  as  a  simple  measure  of 
local  resolution. 

Another  measure  of  system  performance  is  a  measure  of  the  size  of  features  that 
can  be  adequately  resolved  by  the  system.  Chiu  et  al.  (1987)  has  defined  scale  meas- 
urements called  the  horizontal  and  vertical  minimum  resolution  lengths,  //  and  V,  re- 
spectively. Each  represents  the  distance  in  the  respective  dimension  within  which  one 
half  of  the  total  energy  of  the  i":  resolution  kernel  is  confined  under  noise  free  conditions. 
They  are  calculated  from  the  following  expressions, 


/v-i          r,(x",y",z") 
Hx(x',y\  z')  =  2   /2>v2 Et (2>29) 


iiyir'.y. .-')  =  2    2_,*y  - — r (2-30) 


E, 
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'I  2,     v 

r;{x  .  V  ,  r  ) 


n*'.y'.2')  =  2    r^A:2- p (2.31) 


;=i 


where 


(2.32) 


/=i 


is  the  total  energy  of  the  i'h  resolution  kernel  and  r]  is  the  energy  contained  in  the/'1  el- 
ement of  the  /"'  resolution  kernel.  In  three  dimensional  space,  the  center  of  the  box  as- 
sociated with  the  ith  diagonal  element  of  R  is  located  at  (x',y',z'),  and  the  center  of  the 
box  associated  with  the/';  element  of  the  i"'  resolution  kernel  is  located  at  (x",y",  z")  . 
Thus.  A.v.  Ay,  and  Az  are  the  separation  distances  in  each  dimension,  respectively,  be- 
tween these  locations.  A  significant  strength  of  measuring  system  performance  this  way 
is  that  the  RMS  error  and  resolution  do  not  depend  on  the  available  data.  Instead,  error 
and  resolution  can  be  determined  once  the  eigenray  paths  through  the  discretized  ocean 
are  determined  and  the  covariances  of  the  unknown  variable  and  random  error  are 
known. 
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III.     EIGENRAV  PREDICTION  AND  INVERSION  CODE 

DEVELOPMENT 

In  this  chapter,  we  first  discuss  eigenray  path  prediction  for  the  Greenland  Sea  array. 
Eigenrays  are  paths  of  acoustic  signals  extending  between  sources  and  receivers.  Accu- 
rate eigenray  paths  are  needed  to  establish  the  forward  problem.  Errors  in  ray  posi- 
tioning and  travel  time  calculations  introduce  errors  to  the  inverse  solution.  We  then 
discuss  the  modifications  that  have  been  made  to  the  inversion  code  and  their  effects  on 
system  performance. 

A.     EIGENRAV  PATH  PREDICTION 

The  construction  of  the  linear  operator  is  dependent  on  ray  path  information. 
Knowing  the  amount  of  path  length  a  ray  has  in  each  box  is  essential  for  the  estimator 
to  determine  the  distribution  of  unknown  structures.  Collectively,  these  structures  result 
in  a  travel  time  perturbation  datum  as  an  acoustic  signal  passes  from  source  to  receiver. 
We  use  an  algorithm  based  on  ray  theory  to  predict  eigenray  paths  since  ray  theory" 
provides  a  simple  physical  description  and  the  equations  used  in  modeling  are  straight- 
forward. 

The  ray  trace  algorithm  is  a  fourth  order  Runge-Kutta  numerical  integration  tech- 
nique applied  to  Equation  (2.4).  The  depth  of  each  ray  is  calculated  at  1000  range  steps. 
Details  of  the  fourth  order  Runge-Kutta  method  are  given  in  Gerald  (1989).  As  the  ray- 
reaches  significant  points  where  the  vertical  sense  of  motion  reverses,  i.e.  turning  point's 
and  surface  reflections,  the  sign  of  the  second  term  of  Equation  (3.1)  reverses.  (Since 
we  have  not  included  bottom  interacting  eigenrays  in  our  forward  problem,  we  will  not 
consider  bottom  reflections  in  our  discussion  here.)  The  Runge-Kutta  method  cannot 
be  applied  through  a  step  containing  a  significant  point.  Therefore,  in  the  vicinity  of  a 
turning  point  or  a  surface  reflection,  we  apply  a  method  which  assumes  a  locally  con- 
stant sound  speed  gradient  at  depths  near  the  significant  point.  Under  constant  gradient 
conditions,  the  ray  path  is  circular.  From  L"gincius(1970),  we  have  an  expression  for 
curvature  of  the  ray  path. 

x  =  T  [&(v  -  cos  0)  +  gv(2  cos  6v  -  1)]  (3.1) 
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where  gc  is  the  sound  speed  gradient  and  gt  is  the  gradient  of  the  medium  speed.  The 
radius  of  curvature.  0t.  is  just  the  reciprocal  of  k  .  Since  we  use  a  motionless  reference 
ocean  in  the  forward  problem,  the  expression  for  the  radius  of  curvature  simplifies  to 

#■=- C—T  (3-2) 

gc  cos  0 

which  shows  LJt  is  inversely  related  to  the  sound  speed  gradient  with  smaller  radii  re- 
sulting from  stronger  gradients.  In  turning  point  cases,  the  ray  is  projected  through  the 
step  containing  the  point  at  a  constant  ^?,  then  the  Runge-Kutta  integration  continues. 
In  surface  reflection  cases,  the  ray  is  projected  to  the  surface  and  continues  downward 
at  the  same  91,  however,  for  the  downward  path,  the  center  of  the  circle  for  £?  is  sym- 
metrically relocated  about  a  vertical  line  extended  through  the  point  of  reflection.  As  in 
the  turning  point  case,  the  Runge-Kutta  integration  continues  at  the  completion  of  the 
step  until  the  next  significant  point  is  reached.  Figure  2  shows  examples  of  ray  tracing 
through  steps  containing  significant  points. 

To  establish  the  forward  problem,  a  reference  ocean  must  be  chosen.  The  associated 
sound  speed  profile  needs  to  closely  approximate  the  true  conditions  expected  in  the 
Greenland  Sea  so  that  model  predicted  rays  can  be  associated  with  the  data  from  the 
tomographic  array.  When  the  transceiver  moorings  were  deployed,  sound  speed  profiles 
were  generated  at  each  site  from  CTD  data.  The  profiles  indicate  a  thin,  warm  layer  of 
relatively  fresh  water  lies  above  a  sharp  density  gradient.  The  sharp  gradient  extends 
20-50  m  below  the  top  layer.  Below  the  gradient,  the  water  column  approaches  adiabatic 
conditions  at  500-1000  m.  and  is  essentially  adiabatic  below  1000  m.  For  simplicity  in 
modeling,  we  treat  the  sound  speed  profile  as  a  function  of  depth  only.  The  reference 
profile  is  based  on  a  profile  taken  near  the  center  of  the  array.  The  profile  was  taken 
at  the  beginning  of  the  data  collection  period  on  yearday  264  near  Mooring  6  and  is 
shown  in  Figure  3.  Yearday  1  is  defined  as  01  Jan  88  throughout  this  work.  The  pres- 
ence of  the  sharp  and  rapidly  changing  gradients  in  this  profile  near  the  surface  in  con- 
junction with  the  source  and  receiver  locations  proves  to  be  challenging  for  our  ray 
tracing  algorithm. 

Under  certain  circumstances,  small  differences  in  sound  speed  profiles  can  result  in 
important  output  differences  from  our  ray  tracing  algorithm.  The  top  40  meters  of  two 
similar  profiles  with  slightly  different  layer  depths  are  shown  in  Figure  4.  From  40  me- 
ters to  the  bottom  at  3000  meters,  the  profiles  are  identical.  Either  profile  is  represen- 
tative of  the  conditions  near  Mooring  6.    We  have  depicted  the  deeper  layer  (DL)  case 
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Figure  2.       Path  Segment  Geometry  Through  Significant  Points:     a)  Turning  point 
example,  b)  Surface  reflection  example. 
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Figure  3.      CTD  Data  at  Mooring  6  on  Ycarday  264:     From  Worchester  and  Howe 
(1989).  Yenrday  1  =  01/01/8S. 


as  a  possible  result  of  entrainment  of  five  meters  of  the  sharp  gradient  found  in  the 
shallower  layer  (SL)  case  by  mechanical  means. 

Typically,  eigenrays  occur  in  bundles  of  four  distinct  rays  due  to  the  near  surface 
location  of  the  tranceiver  elements  and  near  half  channel  conditions  of  the  Greenland 
Sea.  Two  rays  in  the  bundle  are  launched  downward  within  a  few  tenths  of  a  degree  of 
each  other.  The  deeper  ray  arrives  at  the  receiver  from  below  while  the  shallower  ray 
reflects  off  the  surface  just  prior  to  reaching  the  receiver,  then  arrives  from  above.  The 
other  two  rays  in  the  bundle  are  launched  upward  at  nearly  the  same  magnitude  launch 
angle  as  the  downward  rays.  Again,  the  ray  launched  at  the  steeper  angle  arrives  from 
below  while  the  shallower  ray  reflects  just  before  the  receiver  and  arrives  from  above. 
Since  the  magnitudes  of  the  receiving  angles  are  nearly  the  same  for  all  four  rays  and  the 
bundles  are  several  degrees  apart,  the  bundles  are  identified  by  a  nominal  receiving  an- 
gles. To  identify  ray  paths  as  eigenrays,  our  algorithm  uses  a  tolerance  of  ±  100  meters 
vertical  displacement  from  the  receiver  as  a  threshold.  This  threshold  usually  allows  us 
to  identify  all  four  eigenrays  at  some  nominal  receiving  angle  using  0.1°  separation  in 
launch  angle. 

Eigenray  paths  between  Mooring  4  and  Mooring  5  with  launch  angles  ±6  are 
shown  in  Figure  5  where  7°  <  6  <  14°.  Figure  5a  shows  paths  based  on  the  DL  profile. 
The  SL  case  is  shown  in  Figure  5b.  With  only  a  five  meter  difference  in  layer  depth, 
several  of  the  eigenrays  predicted  for  the  DL  case  are  not  predicted  for  the  SL  case.  The 
rays  of  the  SL  case  reach  the  top  layer  at  lower  angles  and  continue  in  the  layer  for 
longer  ranges  before  returning  to  the  lower  ocean.  The  rays  of  the  DL  case  reach  the 
surface  layer  at  higher  angles  and  have  less  path  in  the  high  speed  water,  which  results 
in  a  significantly  different  path  from  the  SL  case.  Rays  leaving  a  source  at  the  same 
launch  angle.  6  =  —9.03°,  in  the  two  different  environments  are  are  shown  in  Figure  5c. 
In  the  DL  case,  the  ray  is  projected  to  pass  less  than  100  m  from  the  receiver  and  is 
considered  to  be  an  eigenray.  In  the  SL  case,  a  ray  with  the  same  launch  angle  is 
projected  below  the  receiver  by  over  200  meters. 

The  selection  of  predicted  eigenray  paths  used  for  the  inversion  of  data  from  the 
Greenland  Sea  depends  on  which  paths  best  match  the  signal  processing  results  of  the 
array  data.  The  data  we  have  received  for  paths  between  Mooring  4  and  Mooring  5 
indicate,  there  are  four  distinct  paths  at  a  nominal  receiving  angle  of  ±  13°  ,  four  paths 
at  +  10°,  and  the  slowest  arrival  which  occurs  near  0°.  Therefore  we  have  chosen  to 
include  the  eight  steepest  eigenrays  predicted  in  the  DL  case  for  our  inversion  of 
Greenland  Sea  data.    These  rays  are  at  nominal  angles  that  best  match  the  true  ocean 
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Figure  4.      Near  Surface  Sound  Speed  Profiles:     Both  profiles  arc  representative  of 
the  profile  at  Mooring  6. 

data.    In  Chapter  IV,  we  will  continue  our  discussion  of  the  inversion  of  data  from  the 
deployed  array. 

The  presence  of  strong  and  rapidly  changing  gradients  in  our  reference  sound  speed 
profile  causes  another  problem  for  our  ray  tracing  algorithm.  At  relatively  shallow 
launch  angles,  typically  with  0  between  +  5°,  the  algorithm  fails  while  projecting  the  ray 
through  a  turning  point  in  the  strong,  near  surface  gradient.  The  exact  cause  of  failure 
needs  to  be  fully  investigated,  but  a  preliminary  look  indicates  a  higher  degree  of  nu- 
merical accuracy  is  required  and  ray  behavior  in  the  vicinity  of  rapidly  changing  gradi- 
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Figure  5.  Sound  Speed  Effects  on  Eigenray  Paths:  Launch  angles:  7°  <,  0  <,  14° 
a)  DL  case  cigeniays.  b)  SL  case  cigenrays.  c)  Comparison  of  rays 
launched  at  0  =  -  9.03°  using  DL  and  SL  profiles. 


ents  is  not  accurately  depicted.  Based  on  these  results,  we  can  expect  the  accuracy  of 
the  ray  tracing  output  to  degrade  at  shallower  launch  angles,  in  the  presence  of  high 
sound  speed  gradients,  and  in  the  presence  of  rapidly  changing  gradients. 

In  conducting  our  studies  involving  the  entire  array,  we  select  101  eigenrays  along 
the  15  paths  between  moorings  shown  in  Figure  6.  The  geometry  of  the  array  results  in 
three  general  types  of  ray  patterns:  1)  Short  paths  between  outer  transceivers,  2)  long 
paths  between  outer  tranceivers.  and  3)  paths  between  the  outer  tranceivers  and  the 
center  transceiver.   An  example  of  each  is  provided  in  Figure  7. 

More  eigenrays  were  predicted,  but  not  used  for  three  reasons.  First,  we  have  less 
confidence  in  the  accuracy  of  ray  traces  of  shallow  eigenrays.  Also,  initial  results  indi- 
cate that  distinct  travel  time  arrivals  for  shallower  ray  paths  may  not  be  available,  and 
finally,  we  feel  using  these  101  eigenrays  in  our  500  box  discretized  ocean  model  is  suf- 
ficient to  determine  sensitivity  effects  without  taxing  the  limits  of  the  microcomputer 
used  to  conduct  the  experiments. 

B.     VERTICAL  LAYERS  AND  RMS  VARIABILITY  DISTRIBUTION 

The  estimator  arrives  at  a  more  realistic  solution  if  a  priori  information  guides  it 
towards  a  statistically  sensible  solution.  Since  we  generally  expect  to  find  more  strongly 
defined  signatures  of  the  unknown  structures  in  the  upper  portion  of  the  ocean  volume, 
we  wish  to  focus  the  resolving  capability  of  the  estimator  in  this  region.  To  accomplish 
this,  we  make  two  adjustments  to  the  inversion  code.  We  vary  the  thickness  of  the  ver- 
tical layers  rather  than  using  equally  spaced  layers  and  we  describe  the  vertical  variability 
distribution  of  unknown  variable  to  the  estimator  through  the  variable's  covariance 
matrix.  Varying  the  thickness  of  layers  allows  us  to  use  more  boxes  of  smaller  size  to 
describe  the  region  of  greater  variability  without  increasing  the  number  of  discrete  boxes. 
Thus,  the  computational  effort  to  generate  a  solution  is  not  significantly  affected. 

To  show  the  effects  of  our  adjustments  on  model  performance,  we  use  an  ocean 
volume  which  is  240  km  x  240  km  in  the  horizontal,  and  3  km  deep.  The  horizontal  di- 
mensions are  evenly  divided  into  ten  24  km  segments.  There  are  five  vertical  layers. 
Starting  at  the  surface,  the  layers  are  100  m.  200  m,  450  m,  750  m.  and  1500  m  thick, 
respectively.  Because  of  the  uneven  vertical  spacing,  the  lower  boxes  generally  contain 
a  larger  portion  of  eigenray  paths  than  if  the  model  used  evenly  spaced  layers.  As  a  re- 
sult, the  estimator  places  greater  weight  on  the  lower  boxes  when  generating  solutions. 
Figure  8  depicts  the  resolution  kernel  peak  values  along  a  vertical  slice  using  variable 
thickness  layers  (VTL).    Figure  9  shows  the  resolution  kernel  peak  values  for  an  equally 


22 


216   — 


19? 


E 


o 

c 
re 


120 


96 


72    — 


48 


n r 

0      24     4g      72      96     120    144     188    192    2II 

range   (km) 


Figure  6.      Plan  View  of  Array  Geometry 

spaced  (FQS)  discretized  ocean  having  layers  of  600  m  each.  Both  cases  use  the  same 
eigenray  paths  and  the  slices  arc  taken  parallel  to  the  y-axis  along  the  seventh  increment 
in  x,  that  is  for  144km  <  x  <  168km  (sec  Figure  6).  These  slices  essentially  parallel  the 
ray  paths  between  Mooring  1  and  Mooring  3.  Both  estimators  assume  a  constant  RMS 
sound  speed  perturbation  equal  to  2  m/s  throughout  the  ocean. 

The  resolution  kernel  peak  values  arc  a  measure  of  the  resolving  capability  of  the 
estimator  at  a  particular  location.  They  reflect  the  way  the  estimator  views  the  data 
collection  process  as  an  acoustic  signal  travels  along  an  eigenray  path.  In  the  EQS  ap- 
proach, the  estimator  assumes  that  a  greater  portion  of  the  signal  was  accumulated 
above  1000  m  since  most  of  the  eigenray  paths  are  found  in  this  region,  however,  this 
approach  uses  fewer  boxes  to  describe  the  structure  of  the  upper  ocean.  The  VTL  ap- 
proach, in  contrast,  uses  more  boxes  to  describe  the  upper  region  and  thus  gives  more 
detail  on  the  solution  there.  However,  the  V  1'L  approach  also  shifts  resolution  towards 
the  lower  ocean  where  the  system  has  long  paths  passing  through  the  thick  lowest  layer 
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Figure  7.  Eigenray  Paths  along  Vertical  Slices:  a)  Paths  between  Moorings  1  and 
2a.  b)  Paths  between  Moorings  1  and  3.  c)  Paths  between  Moorings  1 
and  6. 
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Figure  8.      Resolution  Kernel  Peak  Values  (VTL  approach):     A  value  of  1.0  implies 
a  box  is  well  resolved. 

of  boxes.  The  increased  detail  on  system  performance  in  the  upper  ocean  and  the  in- 
creased number  of  estimated  points  provided  by  using  smaller  boxes  arc  characteristics 
of  the  VTL  approach  which  we  wish  to  exploit  since  we  expect  greater  variability  in  this 
region.  At  the  same  time,  we  wish  to  limit  the  effect  of  shifting  system  resolution  to- 
wards the  lower  ocean. 

In  the  oceans  described  above,  the  estimators  do  not  expect  greater  variability  in  the 
unknown  field  in  any  particular  region.  'I  his  allows  the  estimator  the  freedom  to  gen- 
erate solutions  which  show  features  of  unusually  high  intensity  in  regions  where  we  re- 
alistically do  not  expect  them.  This  occurs  if  the  solution  meets  the  requirements  of  the 
estimator,  that  is,  it  fits  the  data  and  provides  the  least  mean  square  error.  We  can  in- 
form the  estimator  to  anticipate  greater  RMS  values  of  the  unknown  field  in  the  upper 
ocean  through  the  covariancc  of  the  estimate.  In  this  way  the  estimator  expects  little 
travel  time  perturbation  will  accumulate  in  the  acoustic  wave  as  it  passes  through  the 
deeper  ocean  layers  and  greater  variability  is  encountered  in  the  upper  regions.  As  a 
result,  the  resolution  kernel  peak  values  show  an  upward  shift  in  model  resolution,  off- 
setting the  downward  shift  caused  by  varying  the  thickness  of  layers.  In  Figure  10,  we 
show  a  vertical  slice  of  the  results  of  using  the  VTL  approach  with  a  vertical  dependent 
variability  function.  The  vertical  slice  is  the  same  as  that  used  in  figures  8  and  9.  In 
our  example,  the  RMS  sound  speed  perturbation  is  modeled  as  an  exponentially  decay- 
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Figure  9.      Resolution  Kernel  Peak  Values  (CQS  approach):     A  value  of  1.0  implies 
a  box  is  well  resolved. 

ing  function  with  depth.  An  average  value  is  then  calculated  for  each  layer.  For  the 
example  shown  in  Figure  10,  we  have  characterized  the  function  by  a  surface  RMS  value 
of  6  m's  and  an  e-folding  depth  of  1000  in.  The  average  value  for  the  entire  ocean  using 
this  function  is  1.9  m  s,  which  closely  approximates  the  2  m/s  RMS  5c  used  for  the  ex- 
ample shown  in  Figure  9. 

Though  we  show  an  improvement  in  resolution  in  the  upper  region  of  the  ocean,  it 
does  not  come  without  cost.  As  we  improve  resolution  in  a  region,  we  also  allow  more 
random  noise  to  pass  through  the  estimator.  The  effect  can  be  seen  in  the  standard  de- 
viation of  the  estimate,  and  thus  a  less  reliable  solution.  The  balance  of  higher  resol- 
ution and  lower  random  error  is  a  choice  which  must  be  decided  on  the  basis  of 
resolvable  eigenrays  paths,  array  geometry,  and  the  amount  of  experimental  noise  that 
can  be  expected.  Generally,  we  wish  to  increase  system  resolution  with  a  minimal  in- 
crease in  standard  deviation  of  the  estimate.  Figure  11  shows  the  associated  upward 
shift  in  standard  deviation  of  the  estimate  corresponding  to  an  upward  shift  in  system 
resolution.  The  RMS  Sc  distribution  used  for  the  examples  shown  in  Figures  11a  and 
lib  arc  the  same  as  that  used  for  examples  shown  in  Figures  9  and  10,  respectively. 
Typically,  bias  tends  to  dominate  the  total  error  of  the  system  and  the  increase  in  ran- 
dom error  is  less  important. 
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Figure   10.  Resolution      Kernel      Peak      Values      Using      Depth      Dependent 

Variability:     A  value  of  1.0  implies  a  box  is  well  resolved. 

In  our  modeling  of  the  Greenland  Sea.  we  assume  the  variability  of  sound  speed 
perturbations  and  currents  decrease  exponentially  with  depth,  though  each  variable  may 
be  parameterized  by  a  different  folding  depth.  Reasonable  surface  RMS  dc  values  may 
be  estimated  from  sea  surface  temperature  measurements  derived  from  satellite  data  and 
ship  observations.  Surface  current  variability  may  be  estimated  from  drifting  buoy  and 
CTD  stations  at  the  array-  moorings.  Since  surface  values  may  be  more  reliably  deter- 
mined, the  focus  of  our  sensitivity  study  is  placed  on  the  effect  of  uncertainty  in  sub- 
surface variability.  The  computer  simulated  ocean  on  which  inversions  were  performed 
had  a  surface  RMS  ()c  of  6  m/s  and  a  folding  depth  equal  to  10O0  m.  The  estimator  is 
provided  with  the  same  surface  value  and  then  we  compare  model  performance  using 
folding  depths  of  500  m,  1000  m.  and  3000  m.  A  similar  sensitivity  study  is  conducted 
in  current  tomography.  The  surface  RMS  current  for  the  computer  generated  ocean  and 
the  estimator  is  10  cm  s.  The  ocean  was  generated  using  a  folding  depth  of  700  m. 
System  performance  is  evaluated  using  folding  depths  350  m.  700  m,  and  2100  m  in  the 
estimator.  Figure  12  shows  the  depth  dependent  variability  curves  at  different  folding 
depths.  The  estimator  and  the  simulated  ocean  both  use  an  average  value  based  on 
these  curves  for  each  layer  of  the  discrctized  ocean. 

A  longer  folding  depth  used  in  the  estimator  implies  the  system  expects  greater  var- 
iability in  the  unknown  field  at  depth.    Wc  define  the  half-depth,  %ZS,  as  the  percent 
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Figure  11.  Standard  Deviation  of  the  Estimate  (m/s):  a)  RMS  6c  is  constant 
throughout  ocean  volume,  b)  RMS  5c  is  exponentially  decreasing  with 
depth. 
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Figure   12.       Depth  Dependent  Variability  Curves  at  Various  Folding  Depths:     a) 

Density  tomography  curves,  b)  Current  tomography  curves. 
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of  the  ocean  measured  from  the  surface  which  contains  50%  of  the  integrated  variability 
over  the  ocean  depth.    It  is  expressed  as 

%Z5  =  4—  [ln(0.5  +  0.5e  f  )~|l00%  (3.4) 

where /is  the  folding  depth  and  zec„  =  -3000  m.  is  the  bottom  depth.  Note  that  we  have 
defined /so  that  it  is  positive  in  our  application  here,  that  is/>  0  implies  variability  with 
depth  decreases.  However.  zocn  is  negative  to  be  consistent  with  our  upward  pointing 
positive  z-axis  with  z  =  0  m  at  the  surface.  It  is  interesting  to  note  that  the  half-depth 
depends  only  on  the  folding  depth  and  ocean  depth,  but  not  on  the  surface  value.  This 
allows  us  to  make  comparisons  between  the  two  areas  of  tomography. 

In  Figure  13a.  we  show  the  top  three  layers  of  a  simulated  sound  speed  perturbation 
field  using  the  VTL  approach,  figures  13b,  13c,  and  13d  show  the  corresponding  esti- 
mated fields  based  on  folding  depths  of  500  m,  1000  m,  and  3000  m,  respectively.  The 
simulated  field  has  %Zj  =  21.5%,  or  half  of  the  integrated  variability  is  contained  in 
approximately  the  top  fifth  of  the  ocean.  An  estimator  using  a  folding  depth  at  500  m 
expects  half  of  the  integrated  variability  to  be  contained  in  roughly  the  top  10%  of  the 
ocean.  As  the  acoustic  signal  passes  through  lower  regions  of  the  ocean,  the  signal  en- 
counters higher  variability  than  the  estimator  expects.  The  estimator  interprets  this 
"energetic"'  data  as  the  result  of  more  intense  structures  in  the  upper  ocean  than  are  truly 
present.  The  solution  is  constructed  accordingly.  The  opposite  effect  occurs  when  the 
assumed  folding  depth  is  larger  than  the  true  value.  With  /  =  3000  m.  the  estimator 
assumes  the  ocean  contains  half  the  variability  in  the  top  40%  of  the  ocean  and  tends 
to  smooth  the  variability  over  greater  depths  than  found  in  the  true  solution. 

We  have  made  the  assumption  that  variability  is  an  exponential  distribution.  Under 
these  conditions,  the  estimator  is  more  sensitive  to  underestimating  the  folding  depth 
than  overestimating  it.  Figure  14  shows  the  relationship  between/and  %Zj.  The  slope 
of  the  curve  indicates  system  sensitivity.  For  a  given  folding  depth,  the  slope  is  greater 
when  the  assumed  value  is  less  than  the  true  folding  depth,  indicating  the  system's  in- 
creased sensitivity  to  underestimation.  Since  the  sensitivity  measured  by  %Z5  is  de- 
pendent only  on  folding  depth  and  ocean  depth,  we  can  also  see  from  Figure  14  that 
when  true  folding  depths  are  small,  the  system  is  more  sensitive  to  both  underestimation 
and  overestimation  than  when  the  true  folding  depth  is  large. 

For  current  tomography,  we  find  a  similar  trend.  Figure  15a  shows  a  simulated 
current  field  in  the  top  three  layers  of  the  ocean.    Figures  15b,  15c,  and  15d  show  the 
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Figure   14.       Half-depth  as  a  Function  of  Folding  Depth:     Slope  indicates  model 
sensitivity  to  estimating  true  folding  depth. 

corresponding  estimates  based  on  folding  depths  of  350  m,  700  m,  and  2100  m,  respec- 
tively. As  in  our  density  tomography,  the  estimator  is  more  sensitive  to  underestimating 
the  folding  depth  than  overestimating. 

These  results  of  these  sensitivity  studies  are  based  on  the  assumption  that  the  true 
variability  distribution  is  an  exponentially  decreasing  function  with  depth.  If  the  dis- 
tribution is  significantly  different  from  this,  our  results  may  not  apply,  however,  the 
variation  of  the  half-depth  as  a  function  of  some  other  applicable  parameter  may  still 
be  a  useful  measure  of  sensitivity. 
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Figure  15.  Current  Estimator  Response  Using  Various  Folding  Depths:  a)  Simu- 
lated currents,/  =  700  in.  b)  Estimated  currents,/  =  350  m.  c)  Esti- 
mated currents,/  =  700  m.    d)  Estimated  currents,/  =  2100  m. 
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IV.     PRELIMINARY  INVERSIONS  USING  GREENLAND  SEA  DATA 

The  discussions  from  the  previous  chapters  have  given  us  insight  into  the  expected 
system  performance  in  the  Greenland  Sea  environment.  However,  the  use  of  simulated 
oceans  cannot  explore  even-  aspect  of  system  performance.  Only  using  actual  data  from 
the  tomographic  array  can  disclose  the  extent  of  success  the  system  will  have  as  an  ocean 
monitoring  tool. 

At  present,  a  set  of  travel  time  measurements  from  transmissions  between  Mooring 
4  and  Mooring  5  has  been  processed  for  our  analysis  by  WHOI.  Acoustic  signals  were 
transmitted  six  times  a  day  on  even.'  third  day  starting  in  late  September.  19S8,  and 
ending  in  late  July.  19S9.  covering  the  winter  cooling  season.  The  complete  set  of  data 
for  transmissions  from  Mooring  4  to  Mooring  5  is  used  in  our  preliminary  study.  The 
quality  control  "goodness  estimates"  of  some  data  are.  however,  less  than  optimal.  In 
future  studies,  a  criterion  for  the  minimum  acceptable  goodness  estimate  may  be  set. 
For  our  preliminary  study,  we  use  for  each  ray  an  average  of  all  six  transmissions  to 
determine  the  travel  time  of  the  day.  The  data  indicate  three  bundles  of  eigenrays  were 
resolved  from  signal  processing  at  WHOI.  We  will  make  reference  to  each  bundle  by  the 
magnitude  of  its  nominal  receiving  angle.  The  bundles  consist  of  four  rays  at  13°,  four 
rays  at  10°.  and  one  ray  at  0'.  Using  a  sound  speed  profile  based  on  CTD  taken  near 
the  center  mooring,  our  ray  tracing  algorithm  predicts  four  eigenray  paths'  near  13°and 
four  paths  near  9°  which  we  use  for  our  data  inversion.  Since  the  algorithm  fails  at 
shallow  angles,  the  0°  ray  path  is  constructed  by  an  alternate  method.  We  trace  an  axial 
path  along  the  minimum  sound  speed  axis,  which  is  located  approximately  30  m  above 
the  source  depth.  The  sound  speeds  at  the  source  and  receiver  depths  are  less  than  .3 
m  s  greater  than  the  minimum  sound  speed  of  the  reference  profile.  We  feel  this  is  a 
good  approximation  for  the  0°  ray  since  WHOI  identified  this  arrival  by  using  the  last 
signal  peak  which  was  very  much  bigger  than  noise  (Pawlowicz,1991).  To  display  our 
results,  we  have  taken  a  vertical  slice  in  y-range  as  shown  in  figure  16a.  This  slice  con- 
tains Mooring  4  and  Mooring  5.  The  selected  eigenray  paths,  projected  onto  the  vertical 
slice,  are  shown  in  in  figure  16b. 

To  evaluate  system  performance  and  analyse  the  data,  we  divide  the  ocean  volume 
into  500  boxes  using  the  VTL  layering  scheme  described  in  Chapter  III.  Using  this 
scheme,  the  source  at  Mooring  4  is  in  the  top  ocean  layer  where  48km  <  x  <  72km,  while 
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Figure  16.  Vertical  Slice  and  Eigenray  Paths  for  Moorings  4  and  5:  a)  Orientation 
of  vertical  slice  used  to  display  model  output  and  system  performance 
in  the  vicinity  of  Moorings  4  and  5.  b)  Selected  eigenrays  projected 
onto  the  vertical  slice. 
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the     receiver     at     Mooring     5     is     in     the     second     layer    from    the     top    where 
168km  <  -v  <  192km  We  set  the  a  priori  knowledge  for  the  estimator  as  shown  in  Table 


Table  3.     SYSTEM  PARAMETERS 

Parameter 

Assumed 
Value 

Surface  RMS  Sc 

6  m  s 

Folding  Depth 

1000  m 

Random  Noise 

1  ms 

Correlation  Length  (x.y) 

30  km 

Correlation  Length  (z) 

30o  m 

Once  the  system  parameters  are  established,  the  estimator's  performance  may  be  evalu- 
ated. 

Figure  17  displays  the  resolution  kernel  peak  values.  The  resolution  of  the  system 
is  highest  in  the  vicinity  of  the  moorings  since  this  is  where  the  three  sets  of  rays  con- 
verge. The  higher  resolution  in  the  upper  1000  meters  is  due  to  the  exponential  distrib- 
ution of  variability  and  path  length  contributions  from  both  the  10°  and  13°  rays  in  this 
region. 

Another  measure  of  system  resolution  is  the  minimum  resolution  length.  As  shown 
in  figure  IS.  the  minimum  resolution  length  in  the  y  direction  is  about  60  km  near  the 
top  center  portion  of  the  vertical  slice  and  near  the  source  and  receiver  sites,  and  in- 
creases rapidly  where  there  are  few  or  no  eigenray  paths.  The  relatively  low  values  are 
due  to  the  convergence  of  all  three  ray  groups  in  these  regions.  Ocean  features  smaller 
than  the  minimum  resolution  length  are  not  adequately  resolved  and  will  appear 
smoothed  in  the  estimate.  That  is.  we  expect  the  estimator  to  generate  less  intense  and 
more  spread  out  features  than  the  true  structure  if  the  true  solution  has  a  length  scale 
less  than  the  resolution  length  at  its  location.  Since  we  are  using  a  two  point  system 
which  has  a  north-south  (y-direction)  orientation,  there  is  virtually  no  resolution  of  fea- 
tures in  the  x-direction.  As  more  data  become  available  from  other  moorings,  the  esti- 
mator will  gain  more  information  in  all  dimensions.  Figure  19  shows  the  vertical 
minimum  resolution  length  which  also  tends  to  be  lowest  where  the  three  groups  of  rays 
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Figure   17.       Resolution  Kernel  Peak  Values  for  Greenland  Sea  Data  Inversion:     A 

value  of  1.0  implies  a  box  is  well  resolved. 

converge.  However,  it  is  interesting  to  note  the  unusually  large  minimum  resolution 
length  of  1700  m  in  the  top  center  of  the  slice.  It  appears  to  be  an  anomalous  feature 
caused  by  the  arbitrary  choice  of  the  the  coordinate  system.  The  exact  cause  of  this  lo- 
cally high  value,  however,  needs  further  investigation. 
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Figure   18.       Y-Direction  Minimum  Resolution  Length  (km) 
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Figure   19.      Z-Directiou  Minimum  Resolution  Length  (m) 

The  fact  that  the  system  has  low  resolving  power  is  not  surprising  since  there  are 
only  three  bundles  of  eigenrays  involved.  By  using  depth  dependent  variability  in  the 
estimator,  we  have  concentrated  the  resolution  of  the  system  in  the  upper  ocean.  Since 
RMS  Sc  varies  with  depth,  we  normalize  the  RMS  error  to  show  that  the  system  mini- 
mizes error  in  the  upper  ocean,  the  region  of  interest.  Figure  20  depicts  the  RMS  error 
as  a  percent  of  the  assumed  variability  in  the  layer. 

Producing  sound  speed  perturbation  estimates  with  our  estimator  requires  we  con- 
vert the  travel  time  measurements  from  the  array  to  travel  time  perturbations  with  re- 
spect to  a  reference  ocean.  Once  the  sound  speed  perturbations  have  been  estimated, 
we  can  use  Equation  (2.6)  to  determine  an  estimated  sound  speed  held  of  the  Greenland 
Sea. 

For  each  of  the  nine  eigenrays  identified  by  WIIOl,  a  time  averaged  value  was  cal- 
culated. The  corresponding  travel  time  differences  between  the  nine  time  averaged 
Greenland  Sea  values  and  the  nine  predicted  cigenray  arrivals  based  on  ray  theory  range 
from  200-300  milliseconds.  The  differences  are  largely  due  to  structural  differences  be- 
tween the  model  ocean  used  to  predict  our  ray  paths  and  the  "time  averaged"  Greenland 
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Figure  20.      Percent  RMS  Error 

Sea,  however,  model  errors  in  mooring  positioning  and  in  travel  time  calculations  from 
our  ray  theory  based  algorithm  may  also  contribute  significantly  to  the  differences. 
Based  on  receiving  angles  and  sound  speed  data  taken  at  the  mooring  sites,  we  have  a 
high  degree  of  confidence  in  the  positioning  of  the  ray  paths  within  the  discretized  boxes. 
Even  in  a  worst  case  scenario,  the  200-300  millisecond  differences  would  imply  less  than 
500  meter  differences  in  corresponding  path  lengths  over  total  path  lengths  in  excess  of 
120  km.   These  differences  will  not  affect  the  linear  operator  of  the  estimator  greatly. 

Errors  in  travel  time  calculations  even  as  low  as  tens  of  milliseconds,  however,  may 
contribute  significant  errors  in  the  estimates  since  they  are  treated  by  the  estimator  as 
part  of  the  data.  The  scale  of  S  T  in  the  Greenland  Sea  is  on  the  order  of  tens  of  milli- 
seconds, based  on  differences  between  daily  travel  times  and  the  averaged  travel  time  of 
each  of  the  rays.  Though  we  cannot  eliminate  these  errors,  we  can  separate  the  errors 
from  our  time  series  of  estimates  by  basing  the  travel  time  perturbations  on  the  time 
averaged  ocean.  By  choosing  the  time  averaged  ocean  as  our  reference  state,  we  can 
observe  the  response  of  the  estimator  over  the  winter  cooling  season  without  superim- 
posing the  effects  of  model  travel  time  errors.    The  structure  of  the  time  averaged  ocean 
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is  not  known,  however,  and  may  have  range  dependent  features  of  interest.  We  can  find 
this  structure  with  respect  to  the  known  sound  speed  structure  used  to  predict  our  ray 
paths  by  inverting  the  travel  time  differences  between  corresponding  rays  of  the  model 
ocean  and  the  time  averaged  data  from  the  real  ocean.  From  this  information,  we  can 
find  the  sound  speed  structure  of  the  Greenland  Sea. 

Figure  21  shows  the  estimated  Sc  field  of  the  time  averaged  ocean  referenced  to  the 
ocean  used  to  predict  our  ray  paths.  To  produce  this  field,  we  increased  the  horizontal 
correlation  lengths  of  the  estimator  to  100  km  based  on  the  expectation  that  the  time 
averaged  ocean  is  more  highly  correlated  than  any  of  the  daily  perturbation  fields  we 
estimate.  The  solution  shows  a  reasonable  trend.  The  time  averaged  ocean  is  expected 
to  be  cooler  in  the  upper  ocean  since  the  ocean  used  for  ray  prediction  is  based  on  a 
profile  taken  in  September,  a  relatively  warm  time  in  our  data  collection  period.  How- 
ever, the  accuracy  of  the  predicted  structure  is  affected  by  travel  time  errors  from  the  ray 
tracing  algorithm  and  mooring  positioning  in  the  model.  An  independent  source  of  data 
is  needed  to  verify  the  extent  of  the  effects,  but  positioning  errors  may  be  reduced  as 
more  data  sets  of  the  array  are  included. 

Once  the  system  produces  the  first  estimated  perturbation  field,  there  is  little  in- 
crease in  computational  effort  to  generate  a  time  series  of  solutions  using  the  same  esti- 
mator. To  demonstrate  the  response  of  the  estimator  using  the  Greenland  Sea  data  from 
WHOI,  we  plot  a  time  series  of  the  average  Sc  based  on  the  top  two  model  layers  at  the 
position  of  Mooring  4  and  Mooring  5.  The  time  series  are  shown  in  Figures  22a  and 
23a  respectively.  The  general  "spikiness"  of  the  plots  is  the  result  of  having  a  two  point 
system  which  bases  solutions  on  only  three  groups  of  eigenrays.  It  is  difficult  for  the 
estimator  to  precisely  position  features  with  so  little  information.  The  features  appear 
to  move  more  abruptly  in  sequential  solutions  than  could  be  physically  be  expected,  in- 
dicating instability  of  the  solutions.  The  inclusion  of  data  from  other  moorings  will  in- 
crease the  estimator's  knowledge  of  the  structures  in  the  vicinity  of  Moorings  4  and  5, 
and  the  solutions  are  expected  to  show  greater  stability  over  time.  However,  even  with 
instability  affecting  the  short  term  solution,  there  are  indications  of  trends  at  longer  time 
scales  detected.  The  seasonal  cooling  and  warming  through  the  fall,  winter,  spring  and 
early  summer  time  frame  are  evident  at  both  moorings.  Since  sound  speed  is  directly 
related  to  temperature  (Mackenzie.  19S1),  Sc  decreases  with  winter  cooling  and  increases 
with  spring  and  summer  warming.  Events  at  the  20-90  day  scale  are  extracted  from  the 
time  series  by  apphing  a  time  domain  filter  with  cutoff  frequency  .05  cycles  day  at  each 
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Figure  21.  Sound  Speed  Perturbations  of  the  Time  Averaged  Ocean  (m/s):  The 
sound  speed  structure  of  the  time  averaged  ocean  is  the  sum  of  this  field 
and  the  sound  speed  structure  of  the  ocean  used  to  predicted  eigenray 
paths. 

of  the  point  estimates.    The  time  series  of  the  filtered  solutions  near  Moorings  A  and  5 
are  shown  below  their  respective  unfiltered  versions  in  Figures  22b  and  23b. 

Environmental  temperature  data  taken  at  Moorings  4  and  5  sites  at  a  depth  of  100 
meters  has  been  provided  by  WIIOI  and  are  displayed  in  Figures  22c  and  23c.  Com- 
paring the  filtered  6c  estimates  to  the  respective  temperature  data  provides  evidence  that 
the  two  point  tomographic  system  is  detecting  synoptic  scale  events.  Mooring  5  partic- 
ularly shows  good  correlation  between  temperature  and  the  local  6c  estimates.    Of  ■spe- 
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Sound  Speed  Estimates  and  Temperature  Data  Near  Mooring  4:     a) 

Unfikcred  time  series  of  estimates,  b)  Filtered  time  series  of  estimates 
using  cutoff  frequency  =  .05  cycles/day.  c)  Temperature  data  provided 
bv  W1IOI. 
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Figure  23.  Sound  Speed  Estimates  and  Temperature  Data  Near  Mooring  5:  a) 
Unfihered  time  series  of  estimates,  b)  Filtered  time  series  of  estimates 
using  cutoff  frequency  =  .05  cycles  day.  c)  Temperature  data  provided 
bv  WHO  I. 
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cific  interest  are  the  detection  at  Mooring  5  of  the  strong  cooling  event  near  yearday  360 
and  the  corresponding  responses  to  wanning  and  cooling  from  yearday  390  to  about 
yearday  560. 

As  an  example  of  system  output,  a  sequence  of  seven  be  fields  covering  a  three  week 
period  from  yeardays  462  to  4S0  at  three  day  intervals  is  displayed  in  Figures  24-30.  The 
temperature  record  at  Mooring  5  indicated  that  a  strong  warming  event  followed  by  a 
strong  cooling  event  occurred  during  this  period  while  a  weaker  pair  of  wanning  and 
cooling  events  occurred  near  Mooring  4.  The  output  of  the  estimator  is  consistent  with 
the  temperature  data  and  reflects  these  trends. 
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Figure  24.      Sound  Speed  Perturbation  Estimates  (m/s),  Yearday  462 
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Figure  25.      Sound  Speed  Perturbation  Estimates  (m/s),  Vearday  465 
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Figure  26.      Sound  Speed  Perturbation  Estimates  (m/s),  Yeardny  468 
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Figure  27.      Sound  Speed  Perturbation  Estimates  (m/s),  Yearday  471 
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Figure  28.      Sound  Speed  Perturbation  Estimates  (m/s),  Yearday  474 
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Figure  29.      Sound  Speed  Perturbation  Estimates  (in/s),  Vearday  477 
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Figure  30.      Sound  Speed  Perturbation  Estimates  (m/s),  Yearday  480 
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V.     CONCLUSIONS  AND  RECOMMENDATIONS 

In  this  thesis,  we  discussed  the  advantages  and  problems  of  using  a  ray  theory  based 
algorithm  to  establish  the  forward  problem,  developed  an  inversion  code  which  could 
be  applied  to  Greenland  Sea  data,  and  tested  this  code  using  available  Greenland  Sea 
data.  Here  we  offer  some  conclusions  and  recommendations  for  further  study  in  this 
area. 

A.     CONCLUSIONS 

First,  we  looked  at  the  advantages  and  problems  associated  with  using  an  algorithm 
based  on  ray  theory  to  model  eigenray  paths  for  the  Greenland  Sea  tomography  array. 
Ray  theory  provides  a  simple  depiction  of  eigenrays  paths  and  the  algorithm  is  coded 
with  uncomplicated  mathematics.  However,  CTD  data  from  this  area  shows  an  acous- 
tically difficult  enviroment  may  exist  at  least  during  part  of  the  year  and  degrade  the 
reliability  of  the  predicted  eigenray  paths.  Strong  gradients  and  rapidly  changing  gradi- 
ents in  the  vicinity  of  the  transceiver  elements  can  strongly  affect  model  output.  Rays 
which  are  predicted  to  be  turned  in  the  gradient  or  near  the  sharp  change  in  gradient  are 
strongly  influenced  by  small  changes  in  the  sound  speed  profile  in  this  region.  These 
significant  differences  reduce  the  reliability  of  the  resulting  ray  path.  The  algorithm  of- 
ten fails  when  predicting  paths  for  rays  launched  at  shallow  angles.  Rays  launched  at 
steeper  angles  pass  through  these  regions  at  greater  angles  and  are  less  affected.  These 
ray  paths  are  more  reliable. 

Data  prepared  by  YVHOI  indicates  resolvable  eigenrays  between  Mooring  4  and 
Mooring  5  are  at  large  enough  launch  angles  that  they  are  not  turned  in  the  sharp  gra- 
dient. Thus,  our  algorithm  is  adequate  to  model  the  forward  problem  in  this  environ- 
ment. 

In  our  inversion  code  development,  we  investigated  the  effects  of  using  variable 
thickness  layers  vice  an  equally  spaced  scheme  to  discretize  the  Greenland  Sea  for  the 
estimator.  Using  thinner  layers  in  the  upper  ocean  and  thicker  layers  in  the  lower  ocean, 
we  kept  the  total  number  of  boxes  the  same.  The  advantages  of  using  this  scheme  is  that 
it  provides  more  solution  detail  in  the  region  of  interest  for  the  same  computational  ef- 
fort. However,  since  the  thick  lower  layers  contain  greater  path  length  in  this  layering 
scheme,  system  resolution  shifts  toward  the  lower  ocean.  However,  by  specifying  a 
depth  dependent  variability  in  the  covariance  of  the  unknown  field,  we  tell  the  estimator 
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to  expect  greater  variability  in  the  upper  ocean.  As  a  result  of  this  modification,  system 
resolution  shifts  back  towards  the  upper  ocean.  For  our  purposes,  we  have  modeled  the 
depth  dependent  variability  as  an  exponentially  deceasing  function  with  depth.  If  more 
accurate  information  is  known,  it  could  just  as  easily  be  used. 

Since  the  true  variability  may  not  be  known,  we  investigated  model  sensitivity  to 
uncertainty  in  variability  distribution.  Using  an  exponentially  shaped  variability  dis- 
tribution as  our  basis,  we  looked  at  the  effects  of  underestimating  and  overestimating  the 
characteristic  folding  depth  of  the  distribution.  Defining  the  half-depth  as  the  percent 
of  the  ocean  measured  from  the  surface  which  contains  50%  of  the  integrated  variability, 
we  have  an  objective  means  of  measuring  sensitivity.  The  half-depth  gives  us  an  idea 
of  how  the  estimator  is  distributing  variability  in  the  ocean.  We  have  shown  the  change 
in  half-depth  per  unit  change  in  folding  depth  is  always  greater  when  the  true  folding 
depth  is  underestimated.  We  have  also  shown  from  the  same  relationship  that  at  smaller 
true  folding  depths,  the  model  is  more  sensitive  to  both  underestimation  and  over  esti- 
mation. Half-depths  relationships  can  be  developed  for  other  variability  distributions 
and  applied  in  similar  ways  to  measure  model  sensitivity. 

Our  inversion  method  was  tested  by  using  available  data  from  the  Greenland  Sea 
array  as  input  to  our  model.  The  system  used  nine  eigenray  paths  between  the  source 
at  Mooring  4  and  the  receiver  at  Mooring  5.  The  results  of  our  preliminary  application 
of  the  model  show  evidence  of  detection  of  the  seasonal  cooling  cycle  and  detection  of 
synoptic  scale  events  at  time  scales  of  20-90  days.  The  two  point  system  also  shows 
unrealistically  large  fluctuations  in  the  solutions  at  periods  of  less  than  20  days.  This 
indicates  instability  of  the  solutions,  however,  we  expect  stability  to  increase  as  more 
data  becomes  available  from  other  moorings. 

Our  acoustic  tomography  code  is  ready  for  application  to  Greenland  Sea  data.  The 
estimator  can  provide  concentrated  estimates  in  the  region  of  interest  and  allows  for  a 
depth  dependent  variability  distribution  of  the  unknown  variable.  With  these  modifica- 
tions, our  model  can  use  tomographic  data  to  monitor  changes  in  ocean  structure. 

B.     RECOMMENDATIONS 

CTD  data  has  shown  that  the  near  surface  structure  of  the  Greenland  Sea  is  chal- 
lenging to  model  acoustically.  Though  our  ray  theory  based  algorithm  seems  adequate 
for  our  purposes  here,  the  combination  of  the  near  surface  gradients  and  transceiver  lo- 
cation make  this  algorithm  less  than  ideal.  However,  the  convenience  of  using  a  ray 
theory  based  algorithm  to  model  the  lower  ocean  stills  keeps  it  an  attractive  basis  for 
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sound  propagation  prediction.  We  recommend  that  a  more  accurate  means  of  estab- 
lishing the  forward  problem  in  the  Greenland  Sea  be  investigated,  perhaps  using  a  hybrid 
of  ray  theory  and  normal  modes. 

Based  on  the  promising  results  using  just  a  two  point  system,  continued  analysis  of 
the  Greenland  Sea  using  this  code  and  data  from  the  remaining  14  mooring  paths  of  the 
deployed  array  should  yield  useful  estimates  of  ocean  circulation  and  be  able  to  detect 
larce  scale  convective  events. 
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